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ABSTRACT 

This  thesis  examines  three  methods  for  calculating  the 
100(1-  °^)%   lower  confidence  limits  for  the  reliability  of 
a  K-sized  series  system.   Assuming  that  each  component 
reliability  has  a  Beta  distribution,  identical  posterior 
parameters  A  and  B  are  assigned  for  each  component.   A 
computer  simulation  model  is  then  developed  to  determine 

A 

values  of  R  Tr,\>      The  posterior  parameters  are  then 
s  ,L(«  ) 

A 

converted  to  attributes  data,  and  R  T  ,    .    is  computed 

'       s  ,  L  00 

using  classical  CHI-SQUARE  methods,  and  the  WOODS-BORSTING 
Method.   The  three  values  are  then  compared.   Although 
no  alternative  approximations  are  examined,  the  results 
indicate  that  a  high  degree  of  accuracy  in  computing 

A 

R    t         is  possible  with  the  classical  or  WOODS-BORSTING 
s,L(oO 

methods,  and  that  it  may  not  be  necessary  to  resort  to 

A 

costly  simulation  techniques  to  obtain  values  of  R 

s,L(~) 
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I.   INTRODUCTION 


The  problem  of  determining  lower  confidence  limits  on 
the  reliability  of  a  K-sized  series  system  using  component 
attributes  data  only  has  been  studied  from  the  standpoint 
of  both  classical  methods  and  Bayesian  methods.   One  classical 
method,  the  WOODS-BORSTING  Method,  is  used  in  cases  where 
unequal  sample  sizes  of  each  component  in  the  series  system 
are  tested. 

The  classical  CHI-SQUARE  method  treats  system  reliability 
as  a  constant,  and  computes  the  lower  confidence  limit  by 
solving  the  equation 

A 

R 


^j,2F+2 


=  1  - 


2N 


Ls,LO0 

where  R    ,  t  is  the  100(l-cC)%   lower  confidence  limit,  F 
s,L(*C) 

is  the  total  number  of  failures,  and  N  is  the  common  sample 
size,  with  each  component  having  an  equal  sample  size  for 
testing.   Attributes  data,  obtained  by  mission  testing,  is 
used  to  calculate  F. 

The  classical  WOODS-BORSTING  CW-B)  Method  computes  the 
lower  confidence  limit  for  system  reliability  from  the 
formula 


A 

R   _  f    .  =  exp 

s ,  L  M 


■-H2?]/*i-,M 


Naval  Postgraduate  School  Report,  Special  Project 
No.  SP-114,  A  Method  for  Computing  Lower  Reliability  Limits 
Using  Component  Failure  Data  with  Unequal  Sample  Sizes,  by 
J.  R.  Borsting  and  W.  M.  Woods,  pp.  2-10,  June  1968. 
Hereafter  cited  as  WOODS-BORSTING. 


where  R.  is  the  reliability  of  component  i,  i=l,2,...,K, 

A  K 

and  S  is  a  logarithmic  function  of  "JTr. .   The  term  r  is 
a  function  of  f.  (the  number  of  failure  of  component  i), 
and  n.  (the  sample  size  of  component  i). 

In  the  Bayesian  methodology,  the  reliability  of  each  of 
the  K  system  components  is  considered  to  be  a  random  variable. 
Assuming  that  component  reliability  is  Beta  distributed,  a 
prior  set  of  parameters,  A.  and  B.,  is  assumed.   N  items  of 
each  component  are  mission  tested  (ie.,  tested  for  a  length 
of  time  equal  to  the  required  mission  length),  and  successes 
(s . )  and  failures  (f.)  are  obtained  for  each  component 
i,  i  =  1,2,...,K.   Application  of  these  attributes  data  to 
the  prior  distribution  results  in  the  reliability  of  each 
component  being  Beta  distributed,  with  posterior  parameters 
A.  +  s.,  and  B_  +  f..   Assuming  independence  among  components, 
the  theoretical  system  reliability  is 

r  =TTR1 

The  distribution  of  R   is  then  determined,  and  the  value  of 

s  ' 

A  t-V, 

R    ,  ,  is  the  ea        percentile  of  this  posterior  distribution. 
s,L(eO 

With  systems  becoming  more  costly,  the  Bayesian  approach, 
with  its  desirable  feature  of  a  smaller  required  sample  size, 
is  becoming  more  popular.   Past  mission  test  data  of  tech- 
nically similar  systems  can  be  utilized  by  engineers  and 
scientists  to  arrive  at  a  prior  distribution  of  the  under- 
lying component  reliability.   The  posterior  distribution  is, 
however,  difficult  to  determine.   The  product  of  individual 
reliabilities,  each  Beta  distributed,  presents  formidable 


problems  to  the  statistician,  with  regard  to  type  of  dis- 
tribution, and  mathematical  tractability . 

A  computer  simulation  technique,  using  randomly 
generated  Beta  variables  with  the  desired  posterior  parameter 
values  can  be  used  to  determine  the  100il-oL)%    lower  confidence 
limit  for  system  reliability.   Alternatively,  one  might 
approximate  this  value  by  converting  the  posterior  parameters 
into  "equivalent"  attributes  data,  and  using  the  100(l-«*  )% 
classical  lower  confidence  limit  calculated  directly  from 
the  formulas  given  above  as  the  Bayesian  lower  confidence 
limit.   That  is,  use  known  classical  system  reliability 
estimation  methods  to  approximate  the  final  solution  in 
Bayesian  methodology.   The  accuracy  of  this  approximation 
is  the  basis  for  the  research  detailed  in  this  thesis. 


II.   EXPERIMENTAL  PROCEDURE 

In  order  to  develop  a  source  of  data  with  which 

meaningful  comparisons  could  be  made,  it  was  necessary  to 

make  some  basic  simplifying  assumptions.   It  was  decided 

that  an  assumed  posterior  distribution  for  each  component 

would  be  utilized.   Further,  it  was  decided  that  each 

component  in  the  K-sized  series  system  would  have  identical 

posterior  parameters.   As  a  result,  the  parameters  A.  and 

B   become  simply  A  and  B. 
i 

The  next  step  was  to  develop  a  computer  simulation  which 
could  be  used  to  generate  the  distribution  of  system  relia- 
bility, utilizing  the  given  Beta  parameters  A  and  B.   From 

A 

this  distribution,  the  value  of  R    .  ,  could  then  be 

s,LW 

determined . 

Next,  the  posterior  parameters  would  be.  converted  into 

attributes  data.   The  classical  formula  and  the  WOODS- 

BORSTING  formula  could  then  be  applied  and  corresponding 

values  of  system  reliability  lower  confidence  limits 

determined . 

After  tabulating  these  values  for  R  r  ,    x ,  obtained  by 

&  s  ,  L  (<* )  ' 

means  of  the  three  methods,  according  to  A  and  B  values, 
as  well  as  by  K9  the  system  size,  an  attempt  would  be  made 
to  identify  which  tended  to  give  "more  accurate"  results. 
In  order  to  assist  in  this  identification,  three  other 


values  would  be  calculated.   These  were  R  ,  1-(1-R)K, 
and  F  (total  failures),  where  R  is  the  estimated  mean  of 

A  A 

the  posterior  distribution^  defined  by  R  = 


A+B 


A.   COMPUTER  SIMULATION 

Let  X  be  a  random  variable  with  an  exponential  distribution, 

Let 

B 
-In  R  =  J~X. 

where 

X.,  is  expCA) 

X2  is  exp  (A+l) 


XD  is  exp(A+B-l) 

D 

2 
Then  R  has  a  Beta  distribution,  and 

9 


R  =  exp(-^Xi) 

X*' 


A  computer  program  was  developed  to  generate  a  series  of 

3 
B  random  exponential  numbers.   Using  the  above  theorem  these 

numbers  were  then  summed  and  converted  to  a  single  random 

variable,  with  a  Beta  distribution,  and  parameters  A  and  B. 


2 

Operations  Evaluation  Group,  Center  for  Naval  Analyses 

of  the  University  of  Rochester,  Research  Contribution  No.  79, 
Confidence  Limits  for  System  Reliability,  by  J.  Bram, 
Appendix  A,  15  March  1968. 

3 

Emshoff ,  James  R.  and  Sisson,  Roger  L. ,  Design  and  Use 

of  Computer  simulation  Models,  p.  179,  McMillan,  1970. 


This  procedure  was  then  repeated  K  times,  storing  the  numbers, 
and  multiplying  the  numbers  together  as  they  were  developed. 
The  result  was  a  single  number  for  system  reliability.   The 
complete  procedure  was  then  repeated  1000  times,  storing 
each  value  of  system  reliability.   The  1000  values  were 
then  placed  in  order  from  low  to  high.   This,  then,  was  the 
distribution  of  system  reliability,  based  upon  the  given 

A  +-]_ 

values  for  A  and  B.   R    ,  ,  was  then  the  oL   0il  percentile 

s,L(«*) 

of  this  distribution  (le .  ,  *L    =.10,  then  R(100)  is  the  desired 
estimate).   This  program  was  used  to  generate  values  of 

A 

R      r  r~t\    over   a   range   of   values    for   A,    B,    and   K. 
s,L(**)  to  3       3 

B.   WOODS-BORSTING  FORMULA 

This  formula  provides  a  method  for  computing  lower 
confidence  limits  on  the  reliability  of  a  series  system. 


The  formula  is  repeated  below: 


Vlco  =  e? 


{-sM/*i-4*]} 


where 


A 

s 


K      T 

1 


A'l     L 


*2 

n.       Qi 

2(n,_i)    Qi    +   f^T  2~ 


2nj_-3 


^i   "   nT 


and 


>*> 


2n1-3        A 

Q,    + 

2(n    -1)       ± 


ni 


n.-l 


A  o 

Q- 


2^-3 


A 

Qi 

+ 

Qi 

ni 

n.-l 
i 

^OODS-BORSTING.  pp.  2-8. 


with  a  continuity  correction  to  improve  accuracy. 

The  accuracy  of  the  above  formula  is  a  function  of  the 

5 
amount  of  testing,  relative  to  component  unreliability 

.  F  =  niQi 

A  value  of  F  in  excess  of  5  is  considered  to  be  sufficient 

A 

for  reasonable  accuracy  of  R   _  ,.,  . 

s ,  L  C<) 

C.   CHI-SQUARE  FORMULA 

For  a  series  system  with  independent  components,  the 
system  reliability  is  the  product  of  the  component  re- 
liabilities . 

Rs  ■  JRi 

A 

To  find  the  100(l-°0#  lower  confidence  limit  R   _  ,_.»  for 

s ,  L  {**  J 

R  ,  one  proceeds  as  follows: 

Using  mission  test  data  only  on  component  parts,  with 
equal  sample  sizes,  the  values  of  n.  (sample  size),  f. 
(number  of  failures),  and  s.  (number  of  successes)  for 
component  i,  i=l,2,...,K  are  obtained. 
Then 


and 


K 

s  . 

1 

ni 

ft  . 

s 

■  TTr. 

5 
WOODS-BORSTING,  p.  9. 


10 


Rewriting  gives 


_  f. 


r   =  TTC1-Q>  )  ,    Q±  =  *! 
s  +*  n. 


If  Q.  is  small  this  can  be  written  as 

1  kr  *  , 

A       *_     A  ^~  A 

Rs  "  "  " 


*_   a        v A 

=  TTCi-Q,  )  =  1  -  2.  Q, 


The  term  n.Q.  is  distributed  approximately   PoissonCn. Q. ) 

for  Q.-<.1,  and  n.Q>5.  Therefore 

1         K1  X     K       k  * 

V"n.Q,  =  Tf—PCZn.Q^PCnjFQ.) 

for  equal  sample  sizes  ni~^2= ' ' ' =nk=n* 


If  X  is  distributed  Poisson  (A)  and  there  is  one  observation 

7 

a  2 


on  X,  then 


K 

I 


and  since  £_  f .  is  approximately  Poisson  (n£Q^),  then 


n£Q 


=  %t 


JL-l 

k 


'«*,  2F+2         ,    F    =  21  f. 

U(«*)  2  ^rl       X 


and 


?Qi 


/V  A 

=     *ot,2F  +  2    =    1-R      T  .     . 
u(-)         ^H^—  s,L(oc) 


Therefore 


R                   =    1    _  %^,2F+2 
S»LC>  fe 


Parzen  E.,  Modern  Probability  Theory  and  Its  Applications, 
p.  245,  404,  Wiley,  I960. 

7 
Lloyd,  D.  K.  and  Lipow,  M. ,  Reliability:  Management, 

Methods,  and  Mathematics ,  p.  218,  Prentice-Hall,  1962. 


11 


This  formula  is  accurate  if  the  underlying  assumptions 

used  in  its  derivation  are  valid,  that  is,  Q.  is  "small," 

and  n.Q.  greater  than  5.   A  measure  of  Q   being  small  enough 

can  be  developed  from  the  maximum  likelihood  estimators 

of  R  and  R  . 
i      s 

Recall 


then 


and 

(R)K  =  1  -  K(l-R)      for  identical 

attributes  data 

If  these  two  values  are  in  close  agreement  with  each  other,  then 

A 

one  may  assume  that  R  t  (oc)  >    calculated  from  the  classical 

A 

CHI-SQUARE  formula,  is  an  accurate  estimate.   R,  a  point 
estimate  for  R,  is  the  mean  of  the  Beta  distribution;  that 
is 

A 

R  =   A 


Rs 

K 

K 

i-Z 

A+B 


D.   OBTAINING  ATTRIBUTES  DATA 

In  the  application  of  Bayesian  statistics  with  a  given 
prior  Beta  distribution  of  R,  parameters  A'  and  B'  are 
assumed,  based  on   some  measure  of  how  reliable  a  component 
is.   For  example,  if  no  knowledge  of  the  component  reliability 
is  available,  a  uniform  prior  (A'=B'=1)  is  assigned.   This  tends 
to  be  the  exception  rather  than  the  rule,  since  it  is  equiv- 
alent in  the  classical  sense  of  allowing  only  one  success- 
ful test  before  actual  mission  testing  takes  place. 

12 


In  essence,  the  assignment  of  prior  probabilities  A'  and 
BT  is  equivalent  in  the  classical  sense  of  allowing  A'+B'-l 
tests,  with  A'  successes,  prior  to  actual  mission  testing. 
This  is  true  because  the  oc         percentile  point  of  the  Beta 
distribution  is  the  lower  100  (.1-°*)$  confidence  limit  for  . 
R  when  there  are  A'+B'-l  Bernoulli  trials,  each  with 
probability  of  success  R,  and  A'  of  these  trials  are 

o 

successful.    After  mission  testing  N  items,  with  S  successes 
and  F  failures,  the  posterior  distribution  of  component 
reliability  is  Beta  distributed,  with  parameters  A=A'+S, 
and  B=B'+F.   This  in  turn  is  equivalent  to  having  had 

(A'+S)+CB'+F)-1 

A+B-l     tests, 
with 

A'+S 
=A        successes 
and 

(B'+S)-l 
=B-1        failures 
Thus,  given  posterior  parameters  A  and  B,  the  conversion 
to  attributes  data  is 

Number  of  tests    A+B-l 
Successes         A 
Failures  B-l 


8 
DeGroot,  Morris  H.,  Optimal  Statistical  Decisions,  p.  160, 
McGraw-Hill,  1970.       


13 


III.   PRESENTATION  OF  DATA 

A 

Values  of  R       obtained  by  means  of  the  simulation 
s  ,  L  (w ) 

(SIM)  model,  the  WOODS-BORSTING  (W-B)  Method,  and  the 
CHI-SQUARE  or  classical  CCL)  method  are  presented  on  the 
following  pages.   The  values  of  A,  B,  and  K  used  were: 

A;  50,  100,  150 
B;  .2,   3,   4 
K;  3,  5,  7,  9,  10,  20,  30,  40 
In  addition,  ©6  values  of  .05,  .10  and  .20  were  used. 

To  use  the  tables,  first  find  the  page  with  the  desired 
values  of  A  and  B.   Read  down  the  left  column  until  the 

A 

desired  value  of  K  is  located.   The  values  of  R       are 

s,L(«) 

arranged  to  the  right,  under  the  indicated  ©C  value. 

A 

The  point  estimate,  R,  for  each  component  are  given 
at  the  top  of  each  table.   In  addition,  the  right  side  of 
each  table  lists  values  of  F  (total  failures),  as  well  as 
Rk  and  1-(1-R)K. 

As  an  example  of  how  to  use  the  tables,  suppose  one 

A 

wishes  to  compare  values  of  R   _ ,   „xl  for  A=100,  B=3, 
and  K=7.   Turning  to  the  correct  table  and  reading  down 
the  left  and  across,  one  finds 


A=100    B=3    R=.971 

o*  =.10 
K  .  .  .  SIM   CL   W-B 


7  ...  .76  .80   .81  . 
14 
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IV.   CONCLUSIONS 

A 

In  order  to  present  the  tabular  values  of  R   _  ,..  in  a 

s  ,  L  (•< ) 

more  meaningful  manner,  they  were  plotted  as  a  dependent 
variable,  using  K  as  the  independent  parameter,  for  a  fixed 
value  of  A,  B,  andoi  .   Three  curves  were  plotted,  one  for 
each  method.   In  addition  a  new  variable  A  was  defined  as 


A  = 


Ak 
R   - 


[l-U-R)k] 


R  =   A 


A+B 


and  was  plotted  on  the  same  graph,  depending  on  the  values 
of  A  and  B.   This  was  done  to  show  the  increase  In  the  error 

A 

of  R    ,    .  with  K  due  to  the  simplifying  assumptions  made 
s,LC°<) 

in  developing  the  CHI-SQUARE  formula.   These  curves  are 


shown  in  APPENDIX  A. 

For  the  case  A=50,  B=2,  plots  were  made  for  all  three 
®C  values.   There  was  no  major  difference  observed  in  the 
curves  as  «*  was  varied.   As  expected,  the  curves  were  higher 
for  a  given  value  of  K  as©4  was  increased.   The  plots  were 

A 

looked  at  from  the  point  of  view  of  ranges  of  R,  as  defined 
by  A/A+B.   Two  cases  were  examined,  R  —  .98  and  R  <  .98. 


A.   CASE  I;  R  ^*  .98 

A 

For  values  of  Ri   .98  the  value  of  A  as  a  function  of  K 
increased  gradually.   Both  the  W-B  and  CL  curves  provided 

A 

optimistic  values  of  R   _  ,    . .  using  the  simulation  curve 

s,L («■)>      & 

(SIM)  as  a  basis  for  comparison.   The  classical  curve  gave 
a  smaller  relative  error  except  for  values  of  K  less  than 
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about  10,  where  the  W-B  curve  was  closer  to  the  SIM  curve. 

All  three  curves  showed  the  same  trend  downward  as  K 

increased,  and  all  three  curves  were  reasonably  close  to 
each  other. 

B.  CASE  II;  R  <  .98 

For  values  of  ft  <.98  the  curve  of  R       for  the 

s,L(«<) 

classical  method  crossed  the  SIM  curve,  providing  an  op- 

A 

timistic  prediction  for  R       at  lower  values  of  K  and 

s,L(«*) 

a  pessimistic  value  for  higher  K  values.   This  changeover 
occured  for  values  of  K  as  low  as  5  and  a  high  as  35, 

A 

increasing  as  R  increased.   The  W-B  curve  was  once  again 
consistently  optimistic  compared  with  the  simulation  curve. 
The  A  curve  increased  rapidly  for  R  <"  .98. 

C.  GENERAL  CONCLUSIONS 

Since  the  system  reliabilities  which  would  logically 
be  assumed  to  be  of  importance  to  scientists  and  engineers 
are  those  in  excess  of  .75,  especially  in  view  of  high 
system  cost  considerations,  it  is  concluded  from  the  plots 

A 

that  for  a  value  of  R  L(oO,as  obtained  by  the  simulation, 
in  excess  of  .75,  both  the    classical  and  the  WOODS-BORSTING 
Methods  give  quite  good  approximations  to  the  final  solution 
obtained  in  Bayesian  methodology.   Even  though  both  the 
WOODS-BORSTING  and  the  classical  formulas  yield  optimistic 
values,  the  maximum  error  in  the  range  of  (SIM)  R    ^,  ■£  .75 
is  .10,  for  the  parameters  used. 
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D.   AREAS  OF  FUTURE  STUDY 

The  parameters  used  were  chosen  such  that  each  component, 
when  converted  to  attributes  data,  had  the  same  sample  size 
as  all  of  the  remaining  K-l  components.   This  was  done  to 
ascertain  whether  or  not  the  classical  formula,  valid  only 
for  equal  sample  sizes  among  components,  could  be  used  as 
an  approximation  to  the  final  solution.   For  unequal  sample 
sizes,  the  WOODS-BORSTING  formula  must  be  used  —  the 
classical  formula  does  not  apply.   It  remains  to  be  seen 
whether  or  not  the  estimates  obtained  by  means  of  the 
simulation  model  can  be  accurately  approximated  by  the 
W-B  formula,  with  unequal  sample  sizes,  and  varying  numbers 
of  failures,  for  each  component.   The  K  value  was  not 
examined  beyond  40.   A  trend  is  evident  on  the  curves, 
but  research  should  be  done  for  K>40,  rather  than  a  mere 
extrapolation  of  the  data. 

Parameter  values  studied  for  A  and  B  were  limited  to 
those  presented.   A  posterior  parameter  for  A  in  excess  of 
150  is  necessary  to  yield  high  system  reliability  confidence 
limits  for  values  of  K  in  excess  of  15.   This  area  was 
not  examined. 
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DEFINITION  OF  VARIABLES  IN  COMPUTER  PROGRAM 


VARIABLE 

AA 

BB 

DATA(J) 

DEV 

IX 

K 

SYSREL(I) 

XBAR 


Posterior  parameter  A 

Posterior  parameter  B 

Reliability  of  component  (J) 

Standard  deviation  of  system  reliability 

Random  number  seed 

Number  of  components  in  series 

System  reliability  for  replication  (I) 

Mean  value  of  system  reliability 
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c 
c 
c 
c 
c 


c 
c 
c 


c 
c 
c 

c 
c 
c 


c 
c 
c 

c 
c 
c 


SIMULATION  CF  SERIAL  SYSTEM  RELIABILITY  USING 

BETA  DISTRIBUTED  COMPONENT  RELIABILITIES 

DIMENSION  SYSREL(IOOO) 

DIMENSION  HATA(50) 

DIMENSION    BE150) 

DIMENSION  AA150) 

DATA  AA/50*10D./ 

DATA  BB/50*4./ 

K  =  30 

NB=1000 

AB=NB 

IX=232 

SUM=0.0 

SUMSQ=0.0 

DC  66  I=1,NB 

SYSREL(  I  )  =  1. 

DO  67  J  =  1,K 

A=AA( J  ) 

B=3B(J) 

GENERATE  BETA  DISTRIBUTED  RANDOM  VARIABLES 

CALL  BETA(IX,I Y,A,B,BTA) 
IX=IY 

STORE  VALUES  OF  RANDOM  VARIABLES 

DATA! J) =BTA 

COMPUTE  SYSTEM  RELIABILITY 

SYSREL! I )=SYSREL(I )*DATA(J) 

67  CONTINUE 
SUM=SUM+SYSREL(I ) 
SUMSQ=SUMSQ+SYSREL( I)**2 

66    CONTINUE 

XBAR=SUM/AB 

VAR=(SUMSQ-AB*(XBAR**2) )/(AB-l. ) 
DEV=SQRT(VAR) 

ORDER  THE  REPLICATED  VALUES 

CALL  SHELL  (SYSREL,N3) 

PRINT  OUTPUT 

WRITE! 6, 22 )K 

WRITE (6, 13 )X3AR,VAR,DEV 

SUMB=0. 0 

DO  6  3  L  =  1,K 

SUMB=SUMB+BB(LJ 

68  CONTINUE 
WRITE(6,  14)SU^13 

WRITE (6, 23 )SYSREL(50),SYSREL(100),SYSREL(200) 
FORMAT!  //,T5,'  SUM  OF  B  PARAMETERS*  ,T25 , F8  .2  ) 
FORMAT (• 1» ,15,   'NO.  OF  CCMPONENTS= • , T25 , 1  3 ) 


14 
22 

23 


FORMAT!//, T! 


.05, .10,  .20  PERCENTILE  SIMULATION  REL.* 


1,//,T20,3F7«3) 

13  FORMAT!////, T5, ■ OBSERVED  MEAN ' , T20, F 12. 7, / /, T5, 
l'OBSERVED  VAR. • ,T20,F12. 7,//, T5, • OBSERVED  DEV.« 
2T20,F12.7) 

STOP 

END 
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SUBROUTINE    BET  A  (  I  X  ,  IY  ,  A,  B,  BT  A  ) 
NN=B 
SUM=0.0 

DO    100    11=1, NN 
U=l./A 
IX=IX*65539 
R=0.5+IX*.2328306E-9 
IY=IX 

X=-( ALOG(R) )*J 
SUM=SUM+X 
A=A+1. 
100    CONTINUE 

BTA=1./EXP{ SUM) 

RETJRN 

END 


SUBROUTINE  SH5LL(X,N) 

DIMENSION  X(N) 

M  =  N 
10  M=M/2 

IF(M.EQ.O)  RETURN 

K  =  N-M 

J=l 
20  I  =  J 
30  IM=I+M 

IF(X( I ) .LE.X(IM) )  GO  TO  40 

s=xm 

X(I)=X( IM) 
X(IM)=S 
I=I-M 

IF( I .GE.l)  GO  TO  30 
40  J=J+1 

IF  (J  .LE.  K)  GO  TO  20 

GO  TO  10 

END 
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